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Fixed-energy sandpiles with stochastic update rules are known to exhibit a nonequilibrium phase 
transition from an active phase into infinitely many absorbing states. Examples include the con- 
served Manna model, the conserved lattice gas, and the conserved threshold transfer process. It is 
believed that the transitions in these models belong to an autonomous universality class of nonequi- 
librium phase transitions, the so-called Manna class. Contrarily, the present numerical study of 
selected (l+l)-dimensional models in this class suggests that their critical behavior converges to 
directed percolation after very long time, questioning the existence of an independent Manna class. 
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Self-organized criticality (SOC) was introduced in the 
late 80 's as an attempt to explain the ubiquitous variety 
of scale-invariant phenomena in nature [l L Paradigmatic 
examples are sandpile models (see, e.g., 2|), where sand 
is accumulated on a slow time scale and relaxed in form 
of sudden avalanches on a fast time scale. The interplay 
of slow driving and fast relaxation combined with dissi- 
pation at the boundaries drives such systems towards a 
scale-invariant state without any fine tuning of param- 
eters. As a hallmark of SOC, one observes power-law 
distributed avalanche sizes. 

Later it became clear 0, [j] that SOC is closely re- 
lated to nonequilibrium phase transitions into infinitely 
many absorbing states [5[. To make this relation more 
explicit fixed energy sandpiles (FES) were introduced, 
where grains on a periodic lattice follow the same local 
toppling rules as the original sandpile models without any 
drive or dissipation. One of the best studied models is the 
fixed energy version of the nondetcrministic Manna sand- 
pile model 0, [?|, the discrete conserved Manna model 
(DCMM). The transition in this model was found to be 
different from directed pcrcolation(DP) and other pre- 
viously known universality classes [8j . This point of view 
was bolstered by the discovery that various other mod- 
els with conservation laws such as the conserved lattice 
gas model (CLG), the conserved threshold transfer pro- 
cess (CTTP) and the Maslov-Zhang sandpile 
exhibit the same type of universal critical behavior, con- 
stituting the so-called Manna class (MC) of absorbing 
phase transitions with a conserved field which is consid- 
ered today as firmly established 

The paradigm of an independent Manna class, how- 
ever, has always been overshadowed by several disturbing 
observations [12(. The reported estimates for the criti- 
cal exponents of the MC are quite scattered and show 
anomalous scaling behavior. For example, the decay of 
activity from a homogeneously active state, p a (t) ~ tr a , 
and of the survival probability of an avalanche starting 
with a single seed P s (t) ~ t~ s are characterized by dif- 
ferent critical exponents a ^ 6. On the other hand one 
finds that the order parameter exponent /3 and the sur- 




FIG. 1. 



10 10 10 
t [MCS] 

DCMM: (a) Typical time evolution, where single and 
active sites are marked by black and red pixels, respectively, 
(b) Decay in the supercritical phase for random initial condi- 
tions (dotted lines) and natural initial states (solid lines). 

vival probability exponent (3' are equal within error bars 
fl3j . Moreover, the upper critical dimension d c = 4 and 
the mean-field exponents for the MC 13| are same as 
in DP. In addition, the corresponding SOC models are 
known to be unstable against specific perturbations and 
generically flow to DP 14 1. 

In this Letter, we suggest that an independent Manna 
class does not exist. Instead, we find that the apparent 
nonDP behavior, which was seen in many numerical stud- 
ies, is a transient phenomenon, meaning that all models 
of this type are expected to show a DP critical behavior 
after very long time. 

Discrete Model: To support this hypothesis, we first 
revisit the DCMM, showing that its critical behavior de- 
pends crucially on correlations in the initial state. We 
restrict our study to (1+1) dimensions, where possible 
discrepancies between the MC and DP are expected to 
be most pronounced. The DCMM is defined on a chain 
with L sites and periodic boundary conditions, where 
each site j = 1, 2, . . . , L is occupied by rij particles. Sites 
with rij > 2 particles are declared as active, labeled by 
a flag Sj = 1, while empty and singly occupied sites 
are inactive (sj = 0). The model evolves by parallel 
updates, redistributing particles at active sites indepen- 
dently among randomly selected nearest-neighbor sites 
so that the total number of particles TV = J2j n j i s con ~ 
served [see Fig. Q^a)] ■ The DCMM exhibits a continuous 
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FIG. 2. DCMM: (a) Cumulative background density for ran- . . . 0.1024. Here 4> c = 0.89236 obtained from Fig. Hfb). (b) 
dom initial conditions captured at different simulation times. Corresponding stationary densities, (c) Extrapolation of the 
(b) Decay of the density of active sites for natural initial con- static exponent j3. 



ditions on a system with L = 2 sites. 



transition from an active phase into infinitely many ab- 
sorbing states controlled by the average density of par- 
ticles <fi = N/L in the initial state, where the density of 
active sites p a = (sj) plays the role of an order parameter. 

Random initial conditions: Studying homogeneously 
active initial states with a given density <fi, most authors 
used to distribute iV particles randomly on the chain. 
In the active phase (f> > <f> c one expects the density of 
active sites p a (t) to cross over monotonically from an al- 
gebraic decay to a constant value. Surprisingly, we find 
that p a (t) first reaches a minimum, then increases and fi- 
nally saturates at a stationary value [see Fig. [ljb)] . Such 
a nonmonotonic undershooting is quite unusual for ab- 
sorbing phase transitions. Moreover, the curves cannot 
be collapsed onto a single one, indicating the presence 
of several time scales |l2| . We believe that this circum- 
stance is the origin of various numerical inconsistencies 
reported in the literature. 

Explaining the undershooting: As demonstrated in 
Fig. HJa), the dynamics of active sites takes place on 
a heterogeneous background of immobile particles. For 
random initial conditions (i.e.) this background is highly 
disordered so that the spreading process is expected to 
behave like DP with spatially quenched disorder [lij ]. 
slowing down the decay of p a (t) and lowering its value 
in the active phase. However, the disorder of the back- 
ground is not quenched but gets slowly modified by the 
process itself. We find that this feedback gradually ho- 
mogenizes the background on a very slow time scale, 
leading to a subsequent increase of p a (t)- The gradual 
removal of disorder is visualized in Fig. EJa), where we 
plotted the cumulative sum S(j) = (J2l=i n i) ~ Nj/L 
which measures the excess of particles to the left of site j 
compared to the expected average. As can be seen, the 
pronounced density fluctuations of the random initial 
state (black curve) are gradually leveled out, producing 
an almost flat profile (red line) after very long time. 

Natural homogeneous initial states: In presence of un- 
dershooting, the additional curvature in p a (t) can easily 
lead to an erroneous estimate of 4> c and the critical ex- 



ponents [1 21] . As one of our main results, we show that 
the observed undershooting in the supercritical phase can 
be avoided by preparing natural initial states. Following 
an idea introduced earlier in the context of seed simu- 
lations we first let the process run until it becomes 
stationary after the undershooting. As diffusion is known 
to be the conjugate field in the DCMM [13], we then re- 
activate the system by allowing all the particles to diffuse 
for a single Monte Carlo sweep (MCS), which restores a 
high homogeneous activity without destroying the nat- 
ural long-range correlations in the background. After 
reactivation we measure p a (t) as usual. As shown in Fig. 
[ljb), this resolves the problem of undershooting. 

Figure [5Jd shows that the temporal decay of p a (t) for 
random and natural i.e. is in fact very different. For the 
latter we find the critical point </> c = 0.89236(3), which 
differs significantly from the previously reported estimate 
4> c = 0.89199(5) obtained with random i.e. @,IH- For 
the decay exponent, which was previously estimated by 
a = 0.141(24), we get a much larger value a = 0.159(3) 
for natural i.e. which is in agreement with DP. 

Stationary state: Figure [3^ shows the saturation of 
p a {t) in the supercritical phase using natural initial con- 
ditions. As can be seen in panel (b), the stationary den- 
sity plotted against A = 4> — 4>c displays a slight cur- 
vature. Approaching criticality the effective exponent 
/3 e ff (the slope between adjacent data points) exhibits a 
drift from which we can safely conclude that /? < 0.31. 
Plotting /3 e ff against A b with a heuristically determined 
exponent b = 0.45 in panel (c), the data points seem 
to follow an asymptotically straight line, allowing us to 
extrapolate /3 e s visually to criticality. As indicated by 
the blue arrow, the result is again compatible with DP. 
Likewise, by plotting p a (t)/p a against fA"n and search- 
ing for a data collapse (not shown here) we estimate 
i>\\ = 1.75(5). In addition, we estimated the exponent v±_ 
by measuring the correlation function C(r, A) = {sjSj+ r ) 
in the stationary state, which obeys the scaling form 
C(r, A) ~ pj 1 Q(r/S w± ) near criticality. A data collapse 
in Fig. @|a) gives the estimate v± = 1.095(5). Both ex- 
ponents t-'n and v± are in good agreement with DP. 
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FIG. 4. DCMM : (a) Spatial correlation function in the 
stationary state with system size L = 2 11 . (b) Finite-size 
scaling at criticality with natural homogeneous i.e. for L = 2 8 
to 2 13 and corresponding data collapse. 



Finite-size scaling: The dynamical exponent z = 
v\\/v±_ can be determined by finite-size simulations at 
criticality, where the density of active sites is expected 
to obey the scaling form p a (t, L) = t~ a T (jt) ■ Again a 
conventional data collapse in Fig. HJb) leads to the es- 
timate z — 1.51(5) which is somewhat smaller than the 
corresponding DP value 1.58 but larger than previously 
reported estimates. 

Continuous variant of the model: The discrete dy- 
namics of the DCMM has several shortcomings. On the 
one hand, in finite systems the control parameter is 
quantized in steps of l/L, leading to systematic errors 
which have not been taken into account. On the other 
hand, the background noise is difficult to characterize 
because of its telegraphic nature. For this reason we de- 
fine a different variant of the conserved Manna model, 
where the background is modeled by a continuous vari- 
able. We find that this continuous conserved Manna 
model (CCMM) exhibits a particularly clean DP scaling. 

The CCMM is defined on a one-dimensional lattice 
with L sites and periodic boundary conditions, where we 
associate with each site j a real- valued variable Ej called 
energy. A site is declared as active if Ej > 1. In each 
time step all active sites are synchronously updated by 
redistributing the fractions £j and 1 — £j of the energy 
Ej to the two neighboring sites, where £j <E (0, 1) arc 
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random numbers. Clearly this update rule conserves the 
total energy E = J2j Ej . The model exhibits an absorb- 
ing phase transition when the energy density e = E/L 
crosses a certain threshold value e r - 
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CCMM* 


0.141(24) 
0.159(3) 
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TABLE I. Estimates of the critical exponents in the literature 
and the present work (*) compared to DP. 
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FIG. 5. CCMM: (a) Search for the critical point using natu- 
ral initial conditions (inset) and saturation in the supercritical 
phase in a system with L = 10 sites, (b) Stationary density 
of active sites as a function of A = e — e c . (c-d) Data col- 
lapses for spatial correlations and finite size effects analogous 
to Fig. H 



Numerical results for the CCMM: For randomly dis- 
tributed energies the continuous model shows the same 
undershooting of p a {t) as the discrete model. This can 
be avoided by using natural initial states [see Fig. [5fa)]. 
Searching for the transition point (see inset) we find the 
critical energy e c = 0.65797(1) and the decay exponent 
a = 0.1596(2). Repeating the above analysis in the sta- 
tionary state A = e — e c > 0, we find (5 = 0.277(18) 
[solid line in Fig. EJb)] while farther away from critical- 
ity the local slope increases towards 0.381 (dashed line) 
compatible with earlier results Q. Moreover, the data 
collapse of p a {t)/p a against tA v u (not shown) gives the 
estimate v\\ — 1.74(1), while the data collapse for the 
two-point correlation function in Fig. [5tc) gives an esti- 
mate v± = 1.096(4). A finite-size scaling analysis leads to 
the dynamical exponent z = 1.52(1) [see Fig. GJd)]. Both 
v± and z/|| are in good agreement with the corresponding 
DP values. 

Comparison of scaling functions: In Fig. [S] we overlay 
the data collapses of Figs. [4] and [5] with the appropriately 
shifted scaling functions T and Q of DP (black dashed 
lines) which were determined numerically in a directed 
bond percolation process. As can be seen, the curves 
coincide almost perfectly for both the CCMM and the 
DCMM. A similar coincidence for the off critical scaling 
function can be found in the Supplemental Material [i~2| . 

Seed simulations: In seed simulations, we first let the 
system evolve at a given A < until it reaches a natural 
absorbing configuration [l6| . In the discrete (continu- 
ous) model we then choose a random site j with rij > 
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FIG. 6. Left: Finite-size data collapse for the DCMM and 
the CCMM (shifted) compared with the finite-size scaling 
function of DP (dashed black line). Right: Analogous com- 
parison for the stationary correlation function. 



a consequence, the problem of undershooting can be re- 
solved by using homogeneously active initial states with 
natural background correlations. This leads to a slightly 
different estimate of the critical density and therewith 
to different critical exponents which then turn out to be 
mostly compatible with those of DP. As shown in the 
Supplemental Material [ijj], analogous results were ob- 
tained in other sister models belonging to Manna class, 
namely in the (l+l)-dimcnsional CTTP and the CLG on 
a ladder. All our findings suggest that the Manna class is 
not independent but rather an extension or perturbation 
of DP with nontrivial boundary effects 18 1. 
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FIG. 7. Seed simulations in the CCMM (a) and DCMM (b) 
with L — 10 5 sites. The main figure shows both S and ^ 
versus M together with the best fit over four decades (solid 
lines). The inset shows ^j- ~ (^§) Z , fitted over three decades 
(solid lines). The dashed lines indicate previously measured 
values 



(Ej > 0) and transfer a particle (energy portions Ek/2) 
from randomly chosen site(s) k to the target site j until 
it becomes active. This procedure conserves the number 
of particles (energy) without destroying the background 
correlations. We then measure the mean mass M, the 
mean survival time T and the mean area S of the gen- 
erated clusters. In the stationary state near criticality, 
these quantities are expected to scale as 



M 



T 



S 



(1) 



where 7 = vx. + P\\ — 2/3, t = v\\— /3, and a = v± + V11 — /3. 
Although we could not reproduce DP exponents in these 
scaling laws individually, it turns out that A-indepcndent 
ratios of these quantities do exhibit a clean DP scaling. 
For example, Eq. ([T]) implies the relations 



S ~ M 1+ t , S/T ~ M~ 



ST/M ~ (S/T) 2 . (2) 



As shown in Fig. [3 the exponents /3/"f, and z are es- 
timated by 0.127(7), 0.486(9), 1.591(20) for DCMM, and 
by 0.123(3), 0.476(5), 1.579(50) for CCMM respectively. 
Using the scaling relations 7 = u^ + v^ —2/3 and v\\ = (3/a 
these estimates are consistent with those of Table. fl]and 
in excellent agreement with the corresponding DP values 
/3/7 = 0.1214, i/j_/7 = 0.4816, and z = 1.5807. 

Conclusions: We have shown that the dynamics of 
the conserved Manna model in (1+1) dimensions gradu- 
ally removes the disorder in the background by itself. As 
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In this supplement we provide additional details of the numerical simulations and further argu- 
ments to support our conjecture that the Manna class should be related to DP. We discuss the 
dynamics of the background field, the influence of the initial state, and explain why previous works 
may have underestimated the critical threshold. Moreover, we present further numerical data for 
the conserved threshold transfer process (CTTP) and the conserved lattice gas (CLG) on a ladder. 



MOTIVATION TO DOUBT THE EXISTENCE 
OF AN INDEPENDENT MANNA CLASS. 

In our opinion the following observations indicate 
that the commonly accepted scenario of an independent 
Manna class is probably not correct: 

• All existing evidence for an independent Manna 
class is numerical. 

• While for all other established universality classes 
of absorbing phase transitions the numerical esti- 
mates are stable in two or more digits, the estimates 
for the Manna class are scattered over a wide range. 
For example, the exponent a = 0.141(24) quoted in 
[lj comes with an error of 17%, which is unusual. 

• Seed simulations are notoriously difficult to per- 
form and are plagued by the same initial-state- 
dependence as other models with infinitely many 
absorbing states Q. 

• The Manna class and DP arc known to have the 
same mean field theory and the same upper critical 
dimension d c ~\ [3|-|5|. 

• For homogeneous initial states the density of active 
sites in a critical system is expected to decay as 
Pa(t) ~ t~ a while in seed simulations the survival 
probability decays as P s {t) ~ t~ 5 . In the Manna 
class the exponents a = and 6 = (3' ' /v\\ are re- 
ported to be different. On the other hand, one finds 
j5 = j5' within error bars which is a contradic- 
tion in itself. Contrarily Lee and Lee find that the 
scaling relations hold but then there is a mismatch 
in the exponent v± Q . 

• Some of the critical exponents for the Manna model 
and the CTTP [H] were found to be different. Al- 
though such a splitting of universality may be at- 
tributed to the more deterministic character of the 
updates in the CTTP Q, it is still puzzling that 
only some exponents (/3' are different while 
others (e.g. (3, a, S, z) still coincide in both cases. 



On the other hand, Bonachela and Munoz [8( have con- 
vincingly demonstrated that boundary effects of models 
in the assumed Manna class are distinct from those of 
DP. Any questioning of the Manna class must be com- 
patible with their findings. We will come back to this 
point at the end of this supplement. 



DYNAMICS OF THE BACKGROUND FIELD 

All models, which are believed to belong to the Manna 
class (MC), can be interpreted as a DP-like spreading 
process coupled to a conserved background field (i.e. the 
particle density) , which determines the local rate for off- 
spring production. This interpretation was brought into 
a transparent form by Vespignani et al. [§] , who suggested 
that the essential properties of phase transitions in the 
Manna class are captured on a coarse-grained level by a 
system of two coupled Langevin equations for an activity 
field p a (x,t) and a background field <p(x,t): 

dtPa = rp a - bp 2 a + V 2 p a + °~ \JJ> Jl + Up a (f> , (1) 

d t 4> = DV 2 Pa , (2) 

where r,b,a, D,tu are constants and T)(x,t) denotes an 
uncorrelated Gaussian noise. The first equation is just 
the usual Langevin equation for directed percolation ex- 
cept for the term iop a (j), in which the usual fixed rate for 
offspring production is now replaced by the background 
field <f>(x, t). The second equation describes how the back- 
ground field is reordered in the presence of activity. One 
can easily verify that these equations conserve the total 
number of particles J dx<p(x,t). 

Our work is motivated by the conjecture that the 
spreading process in the first equation evolves superdif- 
fusively with a dynamical exponent z < 2 while the con- 
served background evolves only diffusively, providing ef- 
fectively a disordered background on which the spreading 
process takes place. The spreading process in turn re- 
orders the background, but on a much slower time scale. 
Therefore, depending on the initial state, the disorder in 
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the background may have an influence on the temporal 
evolution over a long time. If we use reactived steady 
state configurations as the initial states, here referred to 
as natural initial conditions, the background disorder is 
assured to be low and the critical behavior of the model is 
that of directed percolation. This suggests that an inde- 
pendent Manna class does not exist, rather fixed energy 
sandpiles may be viewed as DP-processes running on a 
self-flattening background. 



INITIAL STATE DEPENDENCE 



Undershooting 



As demonstrated in Fig. 1 of the Letter, a strongly 
disordered background leads to the phenomenon of un- 
dershooting. Since the dynamics of the background is 
slow compared to the spreading process, the situation can 
be compared with a DP process in presence of spatially 
10h17I]. but with the important dif- 



quenched disorder 
ference that in the present case the disorder is not frozen 
but slowly modified by the process itself. This causes a 
glassy behavior which is known to slow down the decay 
of the density of active sites p a . This makes it plausi- 
ble why the estimates for the exponent a = 0.141(24) 
quoted in [l|, describing the decay p a (t) ~ t~ a , tends to 
underestimate the DP exponent a ~ 0.159 and why the 
reported error bar is so large. 

Although the process gradually homogenizes the back- 
ground, this happens on a much slower time scale so that 
a supercritical system first reaches a quasi-stationary 
state, where the disorder of the background is still present 
on scales above the correlation length. Because of the 
expected non- linear response p a {x) ~ (<f>(%) — 4>c) with 
/3 < 1 this disorder will on average decrease the quasi- 
stationary density of active sites, leading to the observed 
undershooting. However, as time proceeds the disorder 
in the background field slowly disappears so that p a {t) 
increases again and eventually reaches its 'true' asymp- 
totic value. This rebound takes more than two orders of 
magnitude in the simulation time, which - if not properly 
taken into account - may have led to additional errors in 
previous steady-state simulations. 



Explaining undershooting: 
Spreading of localized initial states 
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FIG. 1. Cloggy initial state in the CMM at the critical den- 
sity cj> = <f>c- Left: If all particles are initially concentrated 
in a single hump this region will continue to spread until the 
ends touch each other. Middle: Density of active sites as 
a function of time for various system sizes for cloggy (solid 
lines) and natural initial conditions (dashed lines). Right: 
Unlike natural initial states, which scale with the DP expo- 
nent 2 « 1.58 < 2, the curves for cloggy initial conditions 
collapse if they are plotted against t/L 2 . This confirms that 
the background field spreads diffusively (see text). 



initial states 0. To see this let us consider a "cloggy" 
initial state at the critical density <f> = <f> c , where all par- 
ticles are concentrated in a single hump. For example, 
in the DCMM we could pack N = (j> c L particles in a sin- 
gle dense block by placing rij = 2 particles on all sites 
]={)... N/2 while keeping the rest of the system empty. 
Clearly, the process inside the block is supercritical and 
thus the particles will begin to spread, supporting a high 
level of local activity. The spreading region is expected 
to grow until the two ends meet each other because of 
periodic boundary conditions. At this moment the den- 
sity inside the block, which now covers the whole system, 
becomes critical, and therefore the density of active sites 
p a {t) is expected to break down. As demonstrated in 
Fig. [TJ the typical time scale for this density breakdown 
grows quadratically with the system size L. Therefore, 
this numerical experiment confirms that the reordering 
of the background field, as described by Eq. @, is in 
fact diffusive. 

The implication would be that fixed energy sandpiles 
starting with slightly inhomogeneous initial states at the 
critical point are characterized by two competing length 
scales 



and ~ Vt . 



(3) 



The larger length scale is the usual DP-like corre- 
lation length of the cluster of active sites while ^(t) is 
the scale on which the inhomogeneities have been flat- 
tened out. Above criticality, £i(t) grows until it reaches 
a stationary value £i ~ (</> — 4> c )~' ,± ■ At this moment the 
activity is still concentrated in many small humps, giving 



The phenomenon of undershooting can also be ex- 
plained by demonstrating the diffusive character of the 
background dynamics in the case of highly disordered 



1 We are indebted to referee A for suggesting this explanation. 
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on average a smaller activity. However, the humps con- 
tinue to spread diffusively until they touch each other, 
leading to a subsequent increase of the average activity 
density. The 'true' asymptotic density is reached only 
when all humps have been flattened out. 



Comparing various types of initial states 

Having identified undershooting and related memory 
effects as one of the main reasons for the numerical diffi- 
culties in the conserved Manna model with homogeneous 
initial conditions, it is near at hand to create initial states 
for which the disorder in the background field is low. 

In the Letter we favor natural initial conditions by let- 
ting the system first evolve into a stationary or absorbing 
configuration and then reactivating it by diffusion for a 
single time stepH. There are, of course, several other 
ways to create low-disorder initial configurations. For 
the DCMM we may, for example, create 

• flat absorbing states, where vacant sites are dis- 
tributed dcterministically as equidistant as possible 
by setting = 1 for i = 0, 1 . . . [(f>L\ and zero 
otherwise. This state has an artificially flat profile. 

• random absorbing states^], by first setting n, = 
I with probability <fi and zero otherwise, and then 
adding (deleting) particles at vacant (occupied) 
sites so that the total number of particles is [cf>L\ . 
This state is not flat but compared to random ac- 
tive states the disorder is significantly reduced. 

In both cases the states are then reactivated by diffusion 
for a single Monte- Carlo sweep. 

As shown in Fig. [2J the standard random initial con- 
ditions discussed in the Letter (red curve) lead to a pro- 
nounced undershooting, while the other curves appear 
to be quite similar. However, as can be seen in the in- 
set, reactivated random absorbing states (orange curve) 
still undershoot a little, meaning that the disorder in the 
background field is still strong enough to delay the dy- 
namics. On the other hand, natural and flat initial con- 
ditions do not undershoot. This demonstrates that there 
are many ways to engineer an initial state in such a way 
that undershooting is suppressed. 

We favor natural initial conditions for two reasons. On 
the one hand, there has been a similar debate in the con- 



2 Note that the critical density in the DCMM of about 0.9 is rather 
high so that diffusion for a single Monte Carlo sweep generates 
already a high density of active sites. 

3 We thank referee A for this suggestion. 
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FIG. 2. Comparison of the decay for different initial states 
in DCMM with <p - (f> c = 0.0128. The inset shows a zoom 
of the same data. The expected DP scaling function 1Z (see 
text), shown here as black dashed line, compares well with 
the reactivated natural initial condition. 

text of seed simulations, where it turned out that a nat- 
urally correlated background yields the best numerical 
results 0. On the other hand, the data shown in the 
figure is expected to obey the scaling form 

Pa (t,^) = t- a n(tA^), (4) 

where A = <f> — <f> c is the distance from criticality and 
1Z is a universal scaling function which determines the 
form of the curve in the figure. If, as suggested in the 
Letter, the DCMM model indeed belongs to DP, then 
this scaling function should be the same as in DP. In 
the inset of Fig. [5] we overlaid the scaling function of 
ordinary directed bond DP (black dashed line). As can 
be seen, natural homogeneous initial conditions lead to 
a perfect coincidence. This provides a third example for 
the coincidence with DP scaling functions, in addition to 
the finite-size and two-point correlation scaling functions 
discussed in the Letter. The perfect matching of p a {t) for 
natural initial condition with the scaling function of DP 
also confirms that the reactivation process for a single 
time step does not destroy the natural correlations of the 
stationary state. 



Determination of the critical density 

We suspect that the phenomenon of undershooting 
may have led to a systematic underestimation of criti- 
cal densities in previous works. As an example, we show 
that the commonly accepted critical point of the DCMM 
4>c = 0.89199(5) quoted in Q is incorrect and has to be 
replaced by 0.89236(3) - a value far beyond the previous 
error bars. 
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— (1) random i.e., 


$=0.89236 


— (2) random i.e., 


(^=0.89199 


— (3) natural i.e., 
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DP (this work) 



FIG. 3. Temporal decay of the activity density in the DCMM 
for various types of initial states. The inset shows the raw 
data while the main figure displays the same data divided by 
the expected power law of DP. 

In Fig. [3] we illustrate step by step how it comes to 
such a discrepancy: 

(1) If we would simulate the process with random ini- 
tial conditions at the correct critical point <f> c = 
0.89236(3) the glassy influence of the randomness 
slows down the dynamics and leads to a soft pos- 
itive curvature of the data (black curve in Fig. [3]) 
which looks as if the process was slightly supercrit- 
ical. 

(2) This observation may have tempted previous au- 
thors to 'compensate' this curvature by slightly 
lowering the estimate of <fi c . A good-looking "power 
law" p a ~ t~ a (red curve) is then obtained for 
4> c = 0.89199(5). The corresponding slope (red 
dashed line) is roughly —0.14, compatible with pre- 
viously reported estimates of a. 

(3) However, if we simulate the process with natural 
initial conditions for this value for </> c we would ob- 
tain the blue curve which is clearly subcritical. 

(4) Searching for a straight line with natural initial 
conditions (green curve) one obtains the corrected 
estimate <f> c = 0.89236(3) and a power law decay 
Pa{t) ~ t~ a with the exponent a = 0.159(3) which 
is now compatible with DP. 

To support this point further, we demonstrate in Fig. 0] 
how the data analysis for the stationary state would look 
like if we used instead of the corrected estimate the old 
literature value tj> c = 0.89199(5). Already the raw data 
in panel (b) would exhibit a changing sign in the curva- 
ture, which has not been seen so far in absorbing phase 
transitions. Computing the local slopes in panel (c) one 
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FIG. 4. Estimation of the static exponent j3. This figure illus- 
trates how the data of Fig. 3 in the Letter (black lines) would 
look like when using the literature value <f) c = 0.89199(5) (red 
lines) instead of the corrected estimate cj> c = 0.89236(3). 



can easily see where the old estimate /3 = 0.382(19) in 
Ref. [l| comes from (marked by the red arrow) , but with 
increasing numerical effort the effective exponent turns 
out to grow rapidly as we approach criticality, - a clear 
indication that the estimate of the critical point is incor- 
rect. 



SIMILAR RESULTS IN OTHER MODELS 

Conserved Threshold Transfer Process (CTTP) 

To further back up our numerical results, we inves- 
tigate here the so-called Conserved Threshold Transfer 



Process introduced by Rossi et al. [18( , which is believed 
to belong to the Manna class. The CTTP is defined on 
a hypercubic lattice, where each sites can be empty or 
occupied by one or two particles. Sites with two parti- 
cles are active. The model evolves by random-sequential 
updates, i.e. an active site is randomly selected and the 
two particles are independently moved to a randomly se- 
lected nearest neighbor, provided that there are less than 
two particles at the target site. The CTTP may be in- 
terpreted as a fixed-energy Manna sandpile with a height 
restriction, where the occupation of any site can not ex- 
ceed a value n* = 2. 

While Rossi et al. focused on the two-dimensional case, 
several authors have studied the CTTP in one dimen- 
19r 21 1 , obtaining contradicting results (see Table H} . 



sion 



In particular, all authors report different critical densi- 
ties, indicating that they may have been using different 
variants of the model. Moreover, the critical exponents 
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FIG. 5. CTTP in one dimension. The left figure confirms 
the phenomenon of undershooting for random initial condi- 
tions in the supercritical phase for A = <j> — <j) c = 0.0512 
on a chain with L = 2 16 sites. The right panel demonstrates 
the pronounced difference between random and natural initial 
conditions. The data is divided by the expected DP power 
law. 




obtained in these studies differ significantly. 

In this context it is worthwhile to mention the work by 
R. Dickman in 2006 [22j as a remarkable contribution. In 
our opinion he is the first one to fully realize the numeri- 
cal subtleties in simulations of fixed energy sandpiles and 
their importance in determining the critical point of the 
CTTP. To overcome these difficulties he develops various 
empirical methods to compensate the curvature in the 
data. This leads to exponents which, in his words, "are 
not very different from those of DP" (see Table H| . He 
points out that Kockelkoren and Chate [23| as well as 
Ramasco et al. [H| also obtained numerical results close 
to DP in a related model (e.g. /3 = 0.29(2) compared 
to (3dp = 0.277). Despite this obvious proximity to DP 
he was still convinced that fixed energy sandpiles should 
constitute independent universality class as the value of 
the decay exponent a = @/v\\ would be clearly different. 
Here we argue that this last remaining discrepancy in a 
can be resolved by using natural initial states. 



Redoing parts of the simulations for CTTP [22j us- 
ing natural initial conditions we estimate the critical 
density <fi c = 0.929735(15), which is close to the value 
obtained by Dickman |22j by different correction meth- 
ods. Again we observe the phenomenon of undershoot- 
ing, which disappears when natural initial conditions are 
used (see Fig. [5]). At the critical point p a {t) decays as 
t~ a with a = 0.155(5), in excellent agreement with DP 
values. The saturation values of activity in the supercrit- 
ical regime show a pronounced curvature when plotted 
against 4>~4> c in log scale [Fig. EJa)], as already observed 
by Dickman. A linear fit near the critical point gives an 
estimate of j3 = 0.308(2). This is compatible with simi- 
lar results obtained in [22( using finite size scaling. Both 
values of (3 are not far from /3jjp. 

We complement these studies by a collapse of the 
curves p a (i, A) in the supercritical phase for different 



FIG. 6. Conserved threshold transfer process (CTTP). (a) 
Stationary density of active sites as a function of A = <fi — 
4> c . Inset: Saturation in the supercritical phase in a system 
with L = 10 5 sites, (b) Data collapse according to Eq. ([5]) 
giving v\\ — 1.74(1). (c) Data collapse for spatial correlations 
according to Eq.© for different values of A. (d) Finite-size 
scaling with L — 2 7 ,...,2 12 . The best collapse is obtained 
for z = 1.52(2). 



A = <j> — <p c > according to the scaling form 
P„(t,A) = p ft(iA"ii). 



(5) 



As shown in Fig. pjf b) , the best collapse is obtained for 
v\\ = 1.75(1). Moreover, a data collapse of the spatial 
correlation function C(r, A) for different A according to 
the scaling form 



C(r,A) = p 2 a g(rA- 



(6) 



yields v± = 1.13(1) [see Fig. HJc)]. Finally, a standard 
finite size scaling collapse shown in Fig. |6]Jd) gives an 
estimate of z = 1.52(2). Our estimates of the critical 
exponents are listed in Table U along with those obtained 
in earlier works. The use of natural initial conditions 
leads to a coincidence of the decay exponent a with DP 
and brings v±_ and vn closer to the corresponding DP 
values. 

We note that the mysterious splitting of universality 
observed by Liibeck and Heger stating that the con- 
served Manna model and the CTTP in ID differ in some 
but not all critical exponents, disappears when natural 
initial conditions are used. The same applies to apparent 
violations of scaling. 

The CTTP and the conserved Manna model are very 
similar. In fact, if the height restriction in CTTP is lifted 
by allowing n* — > oo, i.e. by taking n* ~ 20 for all prac- 
tical purposes, the dynamics near the critical density be- 
comes identical with that of the conserved Manna model. 
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Ref. 


4>c 


a 


P 


v \\ 


v±_ 


z 


Liibeck [1,^9] 


0.96929(3) 


0.141 


0.382 


2.452 


1.760 


1.393 


Dickman et. al. [20] 


0.92965(3) 


- 


0.412 


2.41 


1.66 


1.45 


Lee [21] 


0.98285(5) 


0.118 


0.396 


3.36 


2.26 


1.49 


Dickman [22] 


0.92978(2) 


0.141 


0.289 


2.03 


1.36 


1.50 


this work 


0.929735(15) 


0.155(5) 


0.308(2) 


1.74(1) 


1.13(1) 


1.52(2) 


DP 




0.159 


0.277 


1.733 


1.096 


1.580 



TABLE I. Reported estimates for the critical exponents of the one-dimensional CTTP compared with DP. The values have 
been rounded to three decimals. 



It is therefore not surprising that both models exhibit the 
same type of critical behavior. 



Conserved Lattice Gas (CLG) 



from such a site hops stochastically as in DCMM or in 
CTTP. However, the numerical analysis by Oliveira was 
restricted to the estimation of only one exponent, namely 
f3 = 0.40(1). 



Another popular model which is believed to represent 
the Manna class is the so-called Conserved Lattice Gas. 
The CLG is defined on a d-dimensional hypercubic lat- 
tice, where each site is either empty or occupied by a sin- 
gle particle (rij = 0, 1). The particles interact via nearest 
neighbor repulsion. This means that a particle, which has 
at least one neighboring particle, tries to move to one of 
the randomly chosen vacant neighbors. If all neighboring 
sites are occupied the particle is unable to move. Thus, a 
particle is active if it has at least one neighboring particle 
and at least one vacant neighboring site. 

The CLG was originally introduced by Rossi et al. [l8| 
who used a parallel update scheme with an exclusion 
principle. Since this method depends on the sequence 
in which the sites are updated, Liibeck used in a later 
study random-sequential dynamics, i.e. an active site is 
randomly selected and the particle jumps to one of the 
neighboring vacant sites with equal probability. 

In one spatial dimension an active site has at most one 
neighboring vacant site, leading to deterministic hopping 
sequences of active particles. As a consequence the tran- 
sition becomes trivial in this case and is characterized by 
integer exponents [25j . As a way out, Oliveira proposed 
to study the CLG on a ladder |26j. Since sites on a lad- 
der have three neighbors each (see Fig. [7]), some of the 
active sites have two vacant neighbors and the particle 





FIG. 7. Conserved Lattice Gas (CLG) on a ladder. Particles 
with a nearest neighbor are active (red bullets), otherwise 
they are inactive (gray bullets). During the update randomly 
selected active particles move to one of the empty nearest 
neighbors, as indicated by the arrows. 
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FIG. 8. CLG on a ladder : (a) Comparison of decay of 
activity from random and natural initial configurations for 
A = 0.0016, 0.0032 for a system of size L — 10 5 . (b) Deter- 
mination of critical point : p a (t)t a saturates at the critical 
point 4> c = 0.48852(2). The random initial condition shows 
an apparent supercritical feature at the critical point <j> c . 

Following this idea, we numerically analyze the CLG 
on a ladder of linear size L, i.e. 2L sites in total. As 
demonstrated in Fig. Eta), we observe again the phe- 
nomenon of undershooting in the decay of p a (t) start- 
ing from a random configuration. This undershooting 
disappears if natural initial conditions are used. With 
natural initial conditions we obtain the critical point 
<fi c = 0.48852(2), which is a bit higher compared to the 
earlier estimate 0.4755(2) reported in (2o| . 

At the critical point the log scale plot of p a (t) gives the 
slope a = 0.1553(1), which is in excellent agreement with 
DP [sec Fig.|EJb)]. In the supercritical regime, the activ- 
ity density saturates at a steady value p a , which varies 
as (0 - <j) c y with /3 = 0.275(6) [see figures Fig. EJa)- 
(b)]. These saturation curves for different values of 



4 Stochasticity in particle transfer is essential for a non-trivial 
critical behavior; critical behaviour of deterministic models are 
known 12711 to be different. 
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FIG. 9. CLG on a ladder : (a) Saturation of ac- 
tivity in the supercritical regime for A = <f> — (f> c = 
0.00005, 0.0001, 0.0002, . . . 0.0256 in a system of size L = 10 5 . 
(b) Corresponding saturation values plotted against <j) — (f> c ; 
the solid line corresponds to the estimate j3 — 0.275(6). (c) 
Determination of v» from a data collapse following Eq. ([5]) 
using the saturation data for 4> - </> c = 0.0016, 0.0032, 0.0064 
and 0.0128 from (a), giving the estimate i/m = 1.76(2). (d) 
A finite size scaling collapse for L = 1000, 2000, 4000, 8000 
yields z — 1.50(2) and /3/u± — 0.251(4); the inset shows the 
corresponding raw data. 



A = (J) ~ 4>c can be collapsed as usual in accordance 
with Eq. ([5]). This is demonstrated in Fig. M, c ) for 
A = 0.0016,0.0032,0.0064,0.0128 giving the estimate 
v\\ = 1.76(2). Likewise, finite size scaling collapse follow- 
ing p a (t,L) = L-P/ v ^F{t/L z ), shown in Fig.^d), gives 
z = 1.50(2) and fi/v± = 0.251(4). 

A summary of these critical exponents can be found 
in Table UH along with the results of previous studies 
of CLG on a ladder and a modified version of CLG in 



one dimension studied by Fiore and Oliveira 28[, where 
both adjacent active particles may jump simultaneously 
to their respective neighboring vacant sites. 



Ref. 


a 




CLG 1D[25] 


1/4 


14 12 


modified [28] 


0.13(1) 


0.277(3) 2.41 0.223(5) 


ladder [26] 




0.40(1) 


this work 


0.1553(1) 


0.275(6) 1.76(2) 0.251(4) 1.50(2) 


DP 


0.159 


0.277 1.73 0.252 1.58 



CLAIM AND OUTLOOK 

In our work we have shown in the example of the one- 
dimensional conserved Manna model that the commonly 
accepted critical point 4> c = 0.89199(5) is incorrect. This 
means that all critical exponents, which were estimated 
on the basis of this critical point, have to be revisited. 

Using natural initial conditions we find a critical be- 
havior which is the same as or at least very close to DP, 
leading us to the convincing but nevertheless speculative 
conjecture that it actually is DP. This conjecture is based 
solely on numerical results and needs to be substantiated 
further by extensive high-performance simulations on su- 
percomputers. 

Note that so far this DP-conjecture is restricted to FES 
models with periodic boundary conditions in one spatial 
dimension, where the conserved background field flattens 
by itself. Whenever this mechanism is disturbed, e.g. by 
other types of boundary conditions or the presence of ex- 
ternal fields, the background field could play a non-trivial 
role and the resulting behavior will differ from DP. In par- 
ticular the results obtained by Bonachela and Munoz 0] , 
who proposed an efficient method to distinguish DP from 
MC by boundary effects, remain valid. This means that 
FES models are still 'more' than directed percolation in 
the sense that the conserved background leads under cer- 
tain circumstances to additional features which cannot 
be seen in ordinary DP models. However, our results 
suggest that these models do not constitute an indepen- 
dent universality class with an autonomous set of critical 
exponents and scaling functions. 

We note that the observation of a self-flattening back- 
ground field seems to be incompatible with earlier field- 
theoretic findings by Pastor-Satorras and Vespignani (2t| , 
who pointed out that quartic vertices become relevant, 
driving the system away from DP. In this context it would 
be important to study whether the results of the present 
paper, which are restricted to one spatial dimension, can 
also be confirmed in higher dimensions. 

Even if the Manna class turned out to be DP, the sim- 
ple conclusion that "SOC is just DP" would be prema- 
ture. On the one hand, driving and dissipation com- 
municate via the background field, which could play a 
non-trivial role in such situations. On the other hand, 
the relationship between SOC and fixed energy sandpilcs 
as such is still debated in the literature We hope 

that our contribution may stimulate further research in 
this direction. 



TABLE II. Reported estimates for the critical exponents of 
CLG on a ladder are compared with DP. 
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TECHNICAL DETAILS 

The presented simulation results were obtained using a 
64-bit version of the four-tap shift-register random num- 
ber generator introduced by R. M. Ziff [33| as well as ran2 
from the Numerical Recipes in C (34J. In both cases we 
have obtained identical results. 
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